Quantitative Phase Diagrams of Branching and Annihilating Random Walks 
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We demonstrate the full power of nonperturbative renormalisation group methods for nonequilib- 
rium situations by calculating the quantitative phase diagrams of simple branching and annihilating 
random walks and checking these results against careful numerical simulations. Specifically, we 
show, for the 2 A — » 0, A — > 2 A case, that an absorbing phase transition exists in dimensions d = 1 
to 6, and argue that mean field theory is restored not in d — 3, as suggested by previous analyses, 
but only in the limit d — > oo. 



PACS numbers: 05.10.-a, 64.60.Ak, 64.60.Ht, 82.20.- 

Thc non-universal properties of continuous phase tran- 
sitions, both at and out of thermal equilibrium, are much 
more difficult to determine than universal quantities such 
as scaling exponents. The latter are generally accessible, 
even for strongly coupled systems, through perturbative 
calculations near critical dimensions, since renormalisa- 
tion group (RG) transformations can then bring them 
in the vicinity of fixed points corresponding to weakly- 
coupled regimes. On the other hand, non-universal prop- 
erties, such as phase diagrams, depend on the whole RG 
flow, which must be controlled to keep track of all the 
microscopic information. Such calculations are thus gen- 
uinely non-pcrturbative. This is usually tempered by 
the fact that the mean-field (or one-loop) approximation 
seems to capture semi-quantitatively the relevant non- 
universal physics of most equilibrium systems. Hence the 
common wisdom that trying to account for fluctuations 
is needless since they are not expected to qualitatively 
alter phase diagrams. In this Letter, performing non- 
perturbative renormalization group (NPRG) calculations 
checked against numerical simulations, we show that this 
belief is drastically infirmed for a simple and widely stud- 
ied nonequilibrium model, whose universal properties are 
otherwise well-understood. 

We treat the case of branching and annihilating ran- 
dom walks (BARW) where particles of a single species 
A diffuse at rate D and undergo the reactions A 2A, 
2A4 0. Whereas mean-field rate equations predict 
that for any nonzero a, the system always reaches an 
active steady state, early simulations 0, m one and 
two space dimensions have revealed the existence of a 
continuous transition to the empty absorbing state, with 
scaling exponents characteristic of the directed percola- 
tion (DP) universality class Cardy and Tauber |4| 
confirmed, through one-loop RG calculations, that fluc- 
tuations can indeed induce, for any nonzero annihilation 
rate A, a dynamical DP transition in dimensions d < 2. 
Their analysis also suggested that mean field theory is 
recovered for d > 2, implying in particular that a tran- 
sition can no longer occur in d = 3, in agreement with 



numerical results by Takayasu and Tretyakov 0. 

Below we determine the quantitative phase diagrams 

of BARW A 2 A, 2 A — > in dimensions 1 to 6, from 
which it follows that the one-loop phase diagrams are 
qualitatively wrong above d = 2, although perturbative 
RG does lead to the correct universal critical exponents. 
Our approach uses and demonstrates the full power of the 
NPRG method, which allows to calculate — in addition to 
the universal DP-class properties — non-universal quan- 
tities even beyond critical dimensions and weak coupling 
regimes. All our results are confirmed, at the quantita- 
tive level, by numerical simulations carefully tailored to 
approach, in a controlled way, the continuous time limit 
where the analysis is performed. Specifically, we show 
that DP-class transitions occur at any nonzero o~/D in 
all dimensions from 1 to 6, but beyond a minimal thresh- 
old {X/D)th for d > 3. The transition line drifts with 
increasing dimension, and (A/Z?)th seems to grow lin- 
early with d. This suggests that the transition occurs for 
all d and that the mean field picture is realized only in 
infinite dimension. Finally we show that the slope of the 
transition lines can be inferred from the master equation 
of a simple two-site model. 

We first briefly review the field theoretical formula- 
tion of BARW, and the nonequilibrium implementation 
of NPRG methods. The stochastic dynamics of BARW is 
described by a master equation, which governs the time 
evolution of the configurational probability P({rii},t). 
For instance the contribution of the birth and death pro- 
cesses at site i reads: 



dP(rii,i) = a dt[(rii-l)P(ni-l,t) - n, ; P(n,,t)] + 
Ad*[(n i +2)(n i +l)P(n i +2,t)-n i (n i -l)P(n4,t)] (I) 



This master equation can be turned into a path integral 
representation in terms of two fields <f> and 4> H- Fol- 
lowing the generating functional of the correlation 
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functions reads Z\ 



J Vcf>V(f) exp(-S[<f>, 4>]) with 



S[<t>,4>] = J d d rdt[cj)(dt- DV 2 )(j)~ o<$4 

+y/2)^{4)(t> 2 - # 2 ) + \<p 2 4> 2 ] . (2) 

This field theory can be then investigated using the 
NPRG method. This formalism — the effective aver- 
age action method — is a continuous implementation, 
on successive momentum shells, of Wilson's block-spin 
concept. It consists in building a sequence of running 
effective actions in which only fluctuations with mo- 
menta larger than the running momentum scale k are 
averaged. At the scale k = A (A -1 corresponding to the 
underlying lattice spacing), no fluctuation is yet taken 
into account and I\- = a coincides with the microscopic 
action S. At k = 0, all fluctuations are integrated out 
and Tk=o encompasses the macroscopic properties of the 
system as it generates the one particle-irreducible corre- 
lation functions. Tk=o is the analogue of the Gibbs free 
energy T at thermal equilibrium. Tk thus interpolates 
between the microscopic action and the effective action. 
The outcome of this procedure is twofold. First, it con- 
stitutes a RG method in the usual sense, in that the uni- 
versal properties of critical phenomena can be derived 
from the behaviour of the flow in the vicinity of fixed 
points. On the other hand, since the NPRG flow can 
relate microscopic quantities to the large-scale behaviour 
of the system, this procedure embodies a calculation of 
the effective action associated with specific microscopic 
models. It thus enables to compute nonuniversal quan- 
tities, which depend on the microscopic definition of the 
system. In this sense it radically differs from perturba- 
tive RG which loses memory of the microscopic details. 
Moreover, the flow equation of Tfe is nonperturbative, and 
as the approximations used do not rely on the smallness 
of a parameter, they are not confined to weak-coupling 
regimes or to the vicinity of critical dimensions. 

We now briefly sketch the construction of the effective 
action Th that only includes large momentum fluctua- 
tions, and give its nonequilibrium NPRG flow (see 
0,0 for a detailed derivation at thermal equilibrium). 
The low and high momentum fluctuation modes are sep- 
arated by a scale dependent mass- like term ASk (4>, 4 1 ) — 
f q <S>{-q)R k (q)<P T (q) where $ = {(f>,4>), and R k {q) is a 
symmetric matrix with zeros on its diagonal, and a cut- 
off function Rk(q) off its diagonal. This cutoff acts as 
an effective mass Rk(q) ~ k 2 for q <C k that suppresses 
the propagation of the low momentum modes, while not 
altering the high momentum ones, Rk(q) — > for q > k. 
In this work, we use Rk(q) = (k 2 — q 2 )9(k 2 — q 2 ) 0, which 
leads to simple analytical expressions (0(x) denoting the 
step function) . Tk is obtained by a Legendre transform of 
logZfc, modified by an additional contribution ASk(ip, VO 
[tp and if) being the functional derivatives of logZk w.r.t. 
sources] , in order to ensure that the proper limits are re- 



covered at the scales k = and k = A. An exact RG 
equation for the flow of Tk is easily derived and underlies 
the nonequilibrium NPRG formalism Q: 



d fc r fc (v>,vO = ^rr{[: 
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+ R k ]- 1 d k Rk}, (3) 



r (2) 

where T k is the field-dependent matrix of the second 
functional derivatives of Tk, and Tr stands for matrix 
trace over internal indices and integration over the inter- 
nal momentum and frequency. Eq. © is a functional 
equation which obviously cannot be solved exactly. To 
handle it, one has to truncate T^. We here implement a 
standard truncation consisting in the leading order of a 
time and space derivative expansion of Tk, which reads: 

r fc (V,V0= / d d rdt[Z k ^(d t ~D k W 2 )^ + Vk(^,^)}. (4) 

We emphasize that the full functional dependence of 
Vk(tp, ip) is dealt with and this truncation hence embodies 
an accurate description of the steady and uniform con- 
figurations of the system, as supported by allprevious 
studies at thermal equilibrium. We refer to |7l l8 ulfllll| 
for detailed discussions of the approximation schemes. 
The flow equations for Zk, D k and Vk(ip,ip) are given 
in 0. They allow us to determine the phase diagrams 

of BARW A 2A, 2A ^ in any dimension. For 
this, we numerically integrate the flow equations from 
an arbitrary initial ultra-violet scale k — A where Tk=\ 
identifies with the microscopic action (J2J) for given A, a 
and D, and we determine the phase of the system (at 
the final scale k = 0) ensuing from this initial condition. 
Figure^ (lines) shows the resulting phase boundaries for 
d= 1 to 6. 

Before commenting on the phase diagrams, we present 
numerical results which fully confirm them. No existing 
simulations are available for a quantitative comparison. 
Worse, in d = 3, no transition was found In these 
simulations the rates a, A, and D were parameterized by 
a single free variable, a drastic constraint likely to prevent 
them from reaching the absorbing state. Moreover, in 0], 
a strict occupation restriction was enforced ( "fermionic" 
model), at odds with the analytical context. Here, we 
adapt the efficient models introduced recently in |l2( to 
simulate arbitrary BARW processes without strict oc- 
cupation restriction. We first notice that modifying the 

models of P3 to study A A 2A, 2A ±* for d > 3 read- 
ily provides evidence of the existence of DP-class transi- 
tions (not shown). This does not allow, however, a quan- 
titative comparison with our analytical results. Reaching 
this aim requires to resort to a time discretisation At of 
the master equation which reproduces as accurately as 
possible the continuous time evolution. A proper dis- 
cretisation is achieved in the limit of At small enough 
to ensure that AP also remains infinitesimal so that the 
obtained critical rates remain invariant under rescaling 
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At. We proceed as follows. At each time step, all sites 
i undergo a parallel update. The on-site reactions are 
ruled by the three independent rates a, D and A: each of 
the rii particles is tested for branching and for diffusion 
to a nearest-neighbor with respective probabilities a At 
and DAt for each neighbor. Simultaneously, each of the 
ni(rii — 1)/2 possible pairs are candidates for annihilation 
with probability 2AA< following Eq. JJ). The occupation 
number rii is then updated according to the successful tri- 
als and the diffusion moves are implemented. We observe 
that the alteration of critical rates induced by changing 
At to At/10 generally does not exceed a few percent 
as long as the probabilities aAt, 2XAt and DAt remain 
smaller than typically 10 -3 . In this limit, draws such 
that the number of annihilated particles would be larger 
than m are exceedingly rare (and rejected anyway). It 
is only in these rather inefficient-looking regimes that we 
can ensure to approach the continuous time limit. The 
code remains however powerful because only (the typ- 
ically few remaining) non-empty sites are visited. The 
critical rates are determined by tracking the decay of 
the total population starting from fully-occupied initial 
conditions. The critical point is characterized by an alge- 
braic decay with the DP exponent 8 = Pz/v, separating 
saturation (supercritical regime) from (quasi-) exponen- 
tial decay. To obtain the typical accuracy of 2% of the 
results presented in Fig. ^ (symbols), system sizes up to 
2 24 sites, and simulation times of up to 10 8 steps were 
necessary. 

Still, the comparison between the analytical and nu- 
merical phase diagrams requires the prior calibration of 
the axes X/D and a/D. Indeed, as simulations imple- 
ment a discrete lattice version of the master equation, 
the "numerical" rates differ from the analytical ones by a 
spatial continuum limit. We therefore introduce an over- 
all offset parameter C which can be interpreted as the 
ratio Ao/A of the corresponding underlying microscopic 
scales — respectively lattice A and ultra-violet A. The 
ratios X/D and a/D are therefore corrected with regards 
to their scaling dimensions by factors C 2 ~~ d and C 2 (see 
Eq. J2J). The unique offset parameter C is fixed by fitting 
the NPRG thresholds C 2 ~ d (X/ D) th (d) to the numerical 
ones. This produces a very accurate match (Fig. G^a)). 
The resulting rescaled analytical full transition lines then 
also closely match their numerical counterparts on the 
range of rates considered in all dimensions from 1 to 6, 
as displayed in Fig[T] [We have checked that fitting other 
quantities to fix C only mildly changes its value, which 
does not spoiled the agreement between the numerical 
and analytical diagrams.] 

We now comment on the phase diagrams. Their main 
feature is that the phase transition exists in all probed 
dimensions, thus qualitatively invalidating the one- loop 
picture. For d > 3, the transition curves are almost par- 
allel straight lines, crossing the X/D axis at a nonzero 
threshold value (A/D)th- Up to d = 6, this threshold 
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FIG. 1: Phase diagrams of BARW A A 2 A, 2 A A in 
dimensions 1 to 6. Lines present NPRG results, rescaled as 
explained in the text. Symbols follow from numerical simula- 
tions. For each dimension, the active phase lies on the left of 
the transition line, the absorbing phase on the right. 
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FIG. 2: (a): Evolution of the thresholds (A/D) th with the 
dimension; (b) and (c): log plots of the transition line in the 
vicinity of the origin in d = 1 and d = 2 with error bars. 



grows linearly with d at a rate 2.248 ±0.058, as extracted 
from a linear fit (see Fig^a)). Further theoretical calcu- 
lations confirm that (X/ D) t h(d) is nearly linear at least 
up to d — 10. This suggests that (A/D) t h becomes infi- 
nite in the limit d — > oo, so that only the active phase 
remains in this limit. In other words, the mean field 
phase diagram seems to be recovered only at infinite di- 
mension, that is neither in d — 3 (one-loop) nor in d = 4 
(upper critical dimension). Below d = 3, the threshold 
vanishes. The approach of the transition curve to the ori- 
gin is quadratic in d — 1, and exponential in d = 2 with a 
coefficient 11.86 ±0.02 analytically, and 11.67 ±0.15 nu- 
merically (Fig. Hfb) and (c)). This is in close agreement 
with the coefficient An predicted by perturbative RG (see 
and 0). 



4 




5 



1 1 1 ^ 1 ^ 1 

4 5 6 7 8 9 10 

XI D 

FIG. 3: Dashed lines: NPRG phase diagrams in d = 1 to 
d = 5 (from left to right); thick solid line: result from the 
two-site model. 

Wc finally give a simple argument supporting our find- 
ing that in the weak diffusion regime the phase bound- 
aries are almost straight lines. Indeed, an effective transi- 
tion line can be estimated from the master equation itself 
in the limit of large a/D and X/D. If D = exactly, sites 
are decoupled and the time evolution of each site can be 
considered independently. The evolution of the proba- 
bility P(n, t) at a given site is governed in this limit by 
the infinite set of coupled differential equations One 
can readily prove from them that the unique stationary 
solution is {P(n = 0) = 1, P(n > 1) = 0}, which cor- 
responds to the absorbing state. For a/X < 10, we have 
checked, by numerically integrating the master equations 
truncated to n — 100, that this state is the unique attrac- 
tor reached at large times. (Including higher occupation 
number equations does not change this result.) Thus, for 
D = 0, the system always ends up in the absorbing phase 
at least up to a / A < 10 and probably for any finite a/X, 
with a relaxation time r growing with a/X. We now ar- 
gue that the absorbing state remains stable when a small 
diffusion is allowed. Qualitatively, we expect this to hold 
true as long as the diffusion time (2dD)~ 1 (2d standing 
for the number of neighbors), remains much larger than 
r since then the particles on one site die before diffusing. 
To make this statement more quantitative, we consider 
the following model. 

Near criticality, the system is very diluted, so we fo- 
cus on an isolated single particle at site /, and model the 
rest of the lattice by a single, initially empty, neighboring 
site J. We numerically integrate a (30 x 30) truncation 
of the master equation for P({nj,nj},t). Again, the 
system always reaches the absorbing phase, but we can 
assume that an active state can be sustained if the par- 
ticle at site / manages to multiply and spread out to the 
neighbouring site before it dies out. We hence study the 
average occupation numbers fii(t) and nj(t). At given A 



and a, as long as D = 0, nj(t) — 0. When D is increased, 
Tij(t) starts to grow and reaches a maximum at t = t m 
before eventually vanishing. This leads us to set up the 
following criterion: the absorbing state is supposed to be 
destabilized if n.j(t m ) reaches one (while ni(t m ) > 1), 
which defines a critical diffusion rate D c . Probing vari- 
ous a/X yields a nearly-straight transition line with ap- 
proximately the observed constant slope (Fig. |2J). This 
argument, as simple as it is, is by no means a rigorous 
proof, but provides further support to the existence of an 
absorbing state in the weak diffusion regime, in all finite 
dimensions. 

In summary, we have provided a clear and coherent 
picture of the quantitative phase diagrams of BARW 
A ^> 2 A, 2 A ^ for d = 1 to 6. We conclude that 
DP-class transitions probably exist in all space dimen- 
sions, and only occur beyond a minimum decay rate 
(A/D)th for d > 2. Note that it is the very existence of 
this threshold which dooms any perturbative approach. 
The close agreement between our simulations and our 
analytical results firmly roots the ability of the NPRG 
method to probe with great accuracy the nonuniversal 
behaviour of a specific model. This offers a promising 
means of investigation of other reaction-diffusion pro- 
cesses. Different types of BARW are obvious candi- 
dates Q , and future work will be devoted to explore the 
so-called PCPD ("pair-contact process with diffusion") 
2A^r 3A, 2A^ U3. 

B.D. wishes to thank G. Oshanin for fruitful discus- 
sions. The LPTHE is Unite Mixte du CNRS, UMR 7589. 
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